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In this work, we study the spin dynamics of Mni2-acetate molecules in the regime of thermally 
assisted tunneling. In particular, we describe the system in the presence of a strong transverse 
magnetic field. Similar to recent experiments, the relaxation time/rate is found to display a series of 
resonances; their Lorentzian shape is found to stem from the tunneling. The dynamic susceptibility 
is calculated starting from the microscopic Hamiltonian and the resonant structure manifests 
itself also in x{ u )- Similar to recent results reported on another molecular magnet, Feg, we find 
oscillations of the relaxation rate as a function of the transverse magnetic field when the field is 
directed along a hard axis of the molecules. This phenomenon is attributed to the interference of 
the geometrical or Berry phase. We propose susceptibility experiments to be carried out for strong 
transverse magnetic fields to study of these oscillations and for a better resolution of the sharp 
satellite peaks in the relaxation rates. 
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I. INTRODUCTION 

In recent years, numerous experimental results on 
macroscopic samples of molecular magnets, especially 
Mni2-acetate and Fes-triazacyclononane, have drawn at- 
tention to the peculiar resonant structure observed mA he 
hysteresis loops and relaxation time measurements ,L!n3 as 
well as in the dynamic susceptibility £m In this paper, 
we concentrate on Mni2 (shorthand for Mni2-acetate). 
At low temperature, the observed relaxation times r are 
long, up to several months and more, and display a series 
of resonances with faster relaxation as a function of an 
external magnetic field directed along the easy (z) axis of 
the sample. These are considered as signs of macroscopic 
quantum tunneling (MQT) of magnetization. 

Typical experimental samples consist of single crys- 
tals or ensembles of aligned crystallites of identical Mni2- 
molecules. Each molecule has eight Mn 3+ and four Mn 4+ 
ions which, in their ferrimagnetic ground state, have a 
total spin S — 10, see Fig. [j]. Due to strong anisotropy 
along one of the crystalline axes (z-direction), there is a 
high potential barrier 

17(5,) = -ASl - BSf, (1) 

with A/k B ~ 0.54K and B/k B « 0.0011K, between the 
opposite orientations of the spin — ±10); the easy 
axis is the same for all the molecules!] The dipolar inter- 
action between the molecular spins, a possible relaxation 
mechanism, has been found to be weak in Mni2.Ln3 In- 
stead, the observed resonant phenomena are attributed 
to quantum tunneling of single spins - the response be- 
ing magnified by the large number of them - interacting 
with the phonons in the lattice. The role of the hyper- 
fine interactions is still under some controversy and is 
only briefly touched upon in the following. 




FIG. 1. Mni2-acetate molecule. The big circles denote 
manganese and the smaller ones oxygen ions. The arrows cor- 
respond to the orientation of the atomary spins in the S = 10 
ground state. The figure is reproduced from [10] under the 
copyright of EDP Sciences. 

The main features of the experimental findings can be 
understood in terms of two competing relaxation mecha- 
nisms: quantum mechanical tunneling through and ther- 
mal activation over the anisotropy barrier. At high 
temperatures (T > 3K or T > 6K, depending on the 
experimentB) , the spins relax predominantly via 
mal activation due to the phonons in the lattice. E£nL3 
In this regime, the relaxation time follows the Arrhc- 
nius law r = t exp(U / k B T) , where U/k B w 60K de- 
notes the barrier height and tq ss 10~ 8 s the inverse 
attempt frequency.t3 When temperature is lowered to 
2K < T < 3K, the time required by the over-barrier 
relaxation increases exponentially but several of the ex- 
cited states still remain thermally populated. In this 
regime, pairs of states on the opposite sides of the bar- 
rier can be brought to degeneracy by tuning the external 
magnetic field. This enhances the probability to tunnel 
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role of the dipolar 
2 ' and combinations 



across the barrier; the tunneling arises due to crystalline 
anisotropy and possible transverse magnetic field at the 
site of the spin. The tunneling amplitudes are the larger 
the closer to the top of the barrier the states are and, 
consequently, the thermal population of the higher states 
plays a key role in relaxation, cf. e.g. Refs. [KJ-jl^ and 
|l4| . At still lower temperatures, tunneling and the re- 
laxation becomes sensitive to fluctuations in the dipolar 
and hyperfine fields JIj In this paper, we concentrate in 
the regime of thermally-activated tunneling. 

Several authors have investigated the spin dynam- 
ics theoretically with the emphasis ranging from "min- 
imal" models, assuming as simple a spin Hamiltonian 
Hs and a model of the surroundings as .possible (in or- 
der to explain experiments, that is)J13llj'oEjt3 to more 
specific models for investigating 
and/or ky perfipq interactions]! 
of thescEJ'oEjO. The thermally activated relaxation 
has typically been studied using a master equation ap- 
proach ±a d^stxiLcJ-lic. time evolution of the spin density 
matrix.m3ldO'li3'H The susceptibility, on the other 
hand, only been treated within a phenomenological 
model.Be3 

The existing theories have been successfull in explain- 
ing the general features seen in experiments. However, 
several points call for further attention. I) A microscopic 
calculation of the dynamic susceptibility is missing alto- 
gether. 2) The effect of a strong transverse magnetic field 
has not been thoroughly studied and, in particular, not 
in the context of the susceptibility. 3) Several authors 
including ourselves have fourid..a_s£ries of side resonances 
to arise in their calculationsE£lli3a - the fact that these 
peaks are not in general (see Ref. ^5] for exceptions) ob- 
served in experiments is not quite clear. In this work, we 
aim at elucidating these points and present calculations 
for the relaxation rates and susceptibility in a unified lan- 
guage that can be conveniently extended to systems with 
stronger coup ling s. We work with a Hamiltonian similar 
to, e.g., Ref. |f, cf. Eqs. (|), (|), and (|), but intro- 
duce an alternative framework to work with the density 
matrix. It is well known that all the resonances can be 
enhanced by a strong transverse magnetic field but we 
show that the resonances can also be reduced and even 
suppressed: Both the relaxation rate, see also Ref. |2^, 
and susceptibility are found to display significant depen- 
dence on the direction of the transverse field suggesting 
that interference effects of a geometrical phase could be 
observed also in Mni2 and, what is more, do so in the 
regime of thermally-activated tunneling. 

The paper is organized as follows. Part || introduces 
the microscopic model used for Mni2 and a discussion 
on the different interaction mechanisms in the system. 
In part III, we develop a time-dependent description of 



the system in terms of a real-time diagrammatic tech- 
nique. This approach is then used to derive and evaluate 
the kinetic equation and the resulting master equation 
governing the spin dynamics. In part 1IIB, we solve the 



t(H) and the static susceptibility xo(H). T he dy namic 
susceptibility \(uj; H) is calculated in part III C and a 
Kubo-type formula is found. Part IV displays the nu- 



merical results for both t(H) and x( w ! H) accompanied 
with a discussion on the results and their relevance to 
experiments. In part |v| we sum up the work. 



II. SYSTEM 

The spin Hamiltonian of a single Mni2 molecule can 
be written in the form Tis = 7i. z + ■ The first term, 



-AS 2 z -BS 4 z -g^ B H z S z , 



(2) 



with S z being the spin component along the easy axis 
(here the z-direction), describes the part that commutes 
with S z . It consists of the anisotropy terms of Eq. (|l|) 
and a Zeeman term which enables external biasing of 
the energies. The anisotmpy constants have been ex- 
perimentally estimatedorZl as A/k-Q — 0.52 — 0.5£ 



and B/k B = 0.0011 - 0.0013K; the g-factor is 1.9.E3 
The resulting energy levels E m for the eigenstates of 
S z \m) = m\m) together with the potential barrier are 
shown schematically in Fig. E| 
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FIG. 2. Schematic of the energy diagram for 2S + 1 = 21 
eigenstates of TL Z together with the functional form of the 
potential barrier (solid line). The figure also shows examples 
of possible transition processes: spin-phonon coupling may 
change m by ±1 or ±2, while the anisotropy term in TLt 
yields tunnel couplings across the barrier - an example of a 
third-order (in B4) process is depicted in the figure. 



The second term in the Hamiltonian. 



H T = --B 4 (Si 



S-) — 9^B(H X S X + HySy 



(3) 



kinetic equation for the field-dependent relaxation times 



does not commute with S z and gives rise to tunnel- 
ing. The i?4-term arises from crystalline anisotropy, 
B A = (4.3 - 14.4) ■ 10~ 5 K (below we use B 4 = 8.6 • 10" 5 K, 
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but the particular choice is unimportant for the re- 
sults obtained) f5 while the second term is the Zee- 
man term corresponding to a transverse magnetic field 
H± = Hs'md (in spherical coordinates, 9 is the polar 
angle away from the z-axis; the azimuth angle is denoted 
</>: H x = H± cos <f> and H y = H± sin <$>), 

Figure || shows the eigenenergies of TLs as a function of 
the longitudinal magnetic field H z . Away from the res- 
onances, the eigenstates \d) resemble the states |m) and 
are also localized onto the different sides of the barrier 
- the linear field dependence of the eigenenergies stems 
from the Zeeman term in Eq. (|J). Close to the reso- 
nances, H.t couples the to) states across the barrier and 
gives rise to avoided crossings in the energy diagram for 
Ed- The magnitude of these splittings directly gives the 
tunneling strengths. Depending on the states in question 
as well as on the magnitude of B4 (H± = for the mo- 
ment), the splittings are found to vary enormously: from 



states to and m' that are a multiple of four apart are 
coupled. Experimental data do not lend support to such 
a rule, however, but rather suggest that all transitions 
are allowed. It turns out that already a tiny transverse 
field in Eq. (||) is sufficient in achieving this - such a field 
may arise due to dipolar and/or hyperfine interactions 
within the sample, see below, as well as due to uncer- 
tainly in the precise angle between the external field 
and the easy axis of the sample. We assume throughout 
the paper a small constant misalignment angle 9 = 1° 
and H± = |H| • sin 9. In places, we wish to investigate 
the effects of a significantly stronger transverse field and 
state so explicitly. The transverse field increases the tun- 
nel splittings if <fi is close to rar/2 (n is an integer), i.e., 
along the x and y-axes, or leads to oscillations in the split- 
tings if 4> is close to one of the directions 7r(2n + l)/4. 
Below we denote .these special directions as the hard axes 
of the moleculesrJ The tunnel splittings are illustrated 



10 K for the choice B^/k^, = 4.3 ■ 10 K and the states in Fig. |4 as a function of H± and for different angles 



to = ±10 up to almost 2K for B 4 /k B = 14.4 • 10~ 5 K and 
the resonance m = ±2. This upper limit is already of 
the same order of magnitude as the level spacing and in 
fact even exceeds it rendering perturbative calculations 
of the tunneling couplings/strengths somewhat question- 
able. Therefore, even though we first formulate every- 
thing in a general form, independent of the choice of ba- 
sis for the spins, we calculate the actual results working 
in the eigenbasis of Tts- This choice of basis provides the 
additional advantage that it allows us to consider arbi- 
trarily strong magnetic fields. 
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FIG. 3. The eigenenergies of the 21 states as a function 
of H z . The inset shows a blow up of the higher energies at 
low fields - note especially the avoided crossings of the states 
m — m = 4 directly coupled by the i?4-term in Ht- There 
are similar, although smaller, splittings for the lower levels as 
well. The three arrows denote the effective barrier height that 
decreases with increasing H z . 

In absence of H±, there is a selection rule to Ht- only 



see Ref. 
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FIG. 4. The tunnel splittings A m ,_ m as a function of the 
transverse magnetic field; H z = 0. The three types of curves 
correspond to three different angles 0: 0° - dashed, 40° - 
dotted, and 45° - solid. 

The usage of Hs of an isolated spin is based, first of 
all, on the assumption that the molecules indeed reside 
in their S = 10 ground state.E3 This is well justified for 
the experimental temperatures below 3 — 6K - the en- 
ergy required to excite the system to the lowest excited 
state with S — 9 is around 30K. The second assump- 
tion is the absence of interaction. In reality, the spins 
interact with each other via dipolar interaction, with 
the nuclear spins via hyperfine interaction, and with the 
phonons of the surrounding lattice. Experimental evi- 
dence shows that .the dipolar interactions are small for 
T > 2 — 2.4K,liSE3 while the hyperfine interactions pro- 
duce an intriiisip broadening of the order of lOmT to the 
spin statesOEj We n eglect these for the moment and 
return to them in Sec. IV C. 



The spin-phonon interaction is mediated by variations 
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in the local magnetic field induced by lattice vibrations 
and distortions. For tetragonal symmetry, this can be 
generally formulated as, cf. Ref. [lj, 

n S p = gi(e xx - e yy ) <g> (Si - Sy) + -g 2 e X y <8> {S x , S y } 

' 7733 \£xz ®{S X) S z } + e yz ®{S v ,S z }) (4) 

+ 77^4 (U X Z ® {S X ,S Z } + UJy Z (g> {Sy,S Z }) 

where e a p and are the symmetric and antisymmet- 
ric strain tensors, respectively, and gi (i=l,2,3,4) are the 
spin-phonon coupling constants. For these, we adopt the 
values from Ref. [l^: 171 = g^/2 = A and \g 2 \ ~ <?i and 
53 1 ~ 34- I n leading order in gi's, H sp produces transi- 
tions such that for, gi^-terms, Am = ±2 and, for 33,4- 
terms, Am = ±1. The phonons themselves are described 
by 
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(5) 



as a bath of noninteracting bosons. 

To be more quantitative, the phonons are assumed to 
be plane waves with a linear spectrum and three modes - 
two transverse and one longitudinal - denoted by a. The 
elements of the strain tensor, e Q( 3 = (d a up + dffU a )/2 and 
LO a f) EE (d )/2, are defined in terms of the local 

displacement vector 
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(6) 



Here b^V are the bosonic operators for phonons with wave 

ka 

vector k, w? is the correponding frequency, and e^' the 
ath element (of x, y, and z) of the polarization vector; 
N is the number of unit cells and M the mass per unit 
cell. 

Above we considered the energy scales inherent to 
TLs, and the spin-phonon rates, which are found be- 
low (Apps. [X] and ^) to be typically of the order of 
10~ 5 — 10~ 4 K, fall in between the extremes of the tun- 
nel splittings. This value is very small compared to the 
stronger tunneling couplings and seems to suit well for a 
perturbative treatment. On the other hand, for the low- 
lying and weakly coupled states the spin-phonon rates 
may be several orders of magnitude larger than the tun- 
nel splittings and one would expect the tunneling to be 
suppressed. However, the intermediate regime, where the 
tunneling and spin-phonon rates are of the same order of 
magnitude, requires some extra care, see below. 



number of molecules in the samples. Let us define 
M(t) in terms of the reduced density matrix p(t) = 
Tr ph [p tot («)] {p tot (t) is the full density matrix) 

M(t)=gp JB (S z (t)) (7) 
= gpB Tr s [S z p(t)] = g /i B / , m ■ p m , m (t) 

n 1 

The reduced density matrix p(t) describes the spin de- 
grees of freedom in the presence of the phonon reservoir 
- its diagonal elements are just the probabilities for the 
spin to be in the states |m). Our strategy is to start with 
the general formulation 



(S z (t))=Tr[S z (t)-pl° 



(8) 



where p ot = p s p ph denotes the initial density matrix 
of the whole system encompassing the spin and phonon 
degrees of freedom. However, this does not contain the 
interaction between the two: we assume that the cou- 
pling is turned on adiabatically and only enters the time 
evolution of S z (t) (in Heisenberg picture and with the 
convention h = 1) 



(^( t )> = Tr[e + < dt ' W(t ' ) ^e-^ ( *' ) .,to t ] 



(9) 



In order to evaluate Eq. (||) and to describe the dy- 
namics of the system more quantitatively, we introduce 
a real-time diagrammatic language that is applicable to 
any mesoscopic system comprising a part with a finite 
number of states linearly coupled to an external heat (or 
particle) reservoirEj 



A. Diagrammatic language 

Equation (Js|) can be written as a diagram depicted 
in Fig. In this equation, the spin-phonon interaction 
terms can be separated from the noninteracting part of 
the Hamiltonian by shifting to the interaction picture 
with respect to Ho = Hs + Ti-ph (indicated with the sub- 
script I) 

(S z (t))=Tr[p ° i .S z (t)] (10) 

= Tr^r ■fe<o dt " H ^ Sz {t) l Te-< dt ' H ^] 

where (T) T is the (anti-)time-ordering operator and p ot 
has been moved to the front (this is allowed by the in- 
variance of trace under a cyclical permutation of its ar- 
guments) . 



> 



III. DYNAMICS 

The magnetization measured in experiments is the 
molecular magnetization M(t) magnified by the large 



<S z (t)>= p? 
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FIG. 5. Diagrammatic representation of Eq. (|9|), reading 
the expression from right to left and starting at time to at 
the initial state described by po ot , then propagating forward 
(upper propagator) in time to time t, where S z (vertex) is to 
be evaluated, and returning backwards (lower propagator) in 
time to to- In such diagrams, the trace is always implicitly 
assumed. 



FIG. 6. Examples of the diagrams arising after integra- 
tion over phonons. The indices on the forward (backward) 
propagators correspond to the spin states and, in particular, 
the first (second) index of p(t). Note that in general m, m', 
and m" may all be unequal if tunneling is allowed along the 
propagator and if we also consider the nondiagonal elements 
ofp(t). 



We proceed by eliminating the phonon degrees of free- 
dom with the aim to obtain an effective theory for the dy- 
namics of the spin system. We first separate the trace and 
introduce time-ordering Tk along the closed Keldysh- 
contour in Fig. || 

Tr[pl° t -S z (t)} = 

= Tr s {pl Tr phP f T K e- l ^ dt ' n ^ t,)l S z {t)i}- (11) 
The exponential, expanded in powers of 7i S p, reads 

T K e~ l Jk dt W-pI* h = H) m (12) 

rn— 

• / dti / (2*2... / dt m TK_{7i sp (ti)i H sp (<2)l •■■ H s p(i m )l} 
JK JK JK 

with t\ > t-x > ... > t m . Each 7i sp (^)i is repre- 
sented as a vertex on the Keldysh contour. In the 
next step, this expansion together with the explicit ex- 
pression ([|) for Ti sp is inserted into Eq. ([ll]) and the 
trace over the phonons is performed by using Wick's 
theorem. As a consequence, the phonon operators are 
pairwise contracted: Tip^p^b^ (t)b(t') = (&t(f)6(f')) and 
Tr ph pP h 6(t)6 t (t') = (b(t)tf(t')). Here, the phonons are 
considered as a reservoir that is not perturbed by the 
single spin - the bath remains in equilibrium and the 
contractions are given by the Bose-Einstein distribution 
(bt 6- ) = [exp(/3cJg CT ) — l] -1 . In terms of the diagrams 
this means that the vertices get coupled (in all possible 
ways) by interaction or phonon lines which correspond 
to just these contractions, cf. Fig. ^. The general rules 
for evaluating the diagrams are given in App. ^ and the 
contributions from the interaction lines are calculated in 
App. § 
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All the diagrams, e.g., those in Fig. |(], are composed 
of two kinds of elements that can be distinguished by 
drawing a vertical line through the diagram and checking 
whether it cuts phonon lines or not. If it does not, the 
corresponding part of the diagram represents free evo- 
lution of the spin; if it does, the part of the diagram 
between successive periods of free evolution corresponds 
to spin-phonon interaction and, more particularly, to a 
transition rate between different spin states. The sum 
of all the possible diagrams with the interaction lines is 
denoted as S and, when evaluated, its terms are 0(g 2n ) 
where g is the spin-phonon coupling constant and n the 
number of interaction lines in the diagram, cf. Appen- 
dices ^ and [b|, and below. The full time evolution of 
(S z (t)) can be expressed as in Fig. ^|a in terms of the two 
kinds of contributions, just discussed. In lowest order in 
g, £ is given by the diagrams shown in part (b) of the 
figure. 
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FIG. 7. Diagrammatic expression for (a) the expectation 
value of S z in terms of the element E, which corresponds to 
the phonon-induced transition rate of the spin, (b) Diagram- 
matic representation for E in 0(g 2 ). 



B. Kinetic equation 



The reduced density matrix defined as 
p m , m ,(t) = ((\m'){m\){t)) 



(13) 



?q m' m'+l m+1 m t 



accounts for the full time evolution of the spin in the pres- 
ence of the spin-phonon interaction. Its diagrammatic 
expansion is analog to Fig. 0a with S z being replaced by 
\m')(m\. Setting t — and differentiating with respect 
to time, we obtain the kinetic equation 
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m + i[Hs,p(t)]= f dt'X(t,t')p(t') (14) 
Jo 

The commutator on the l.h.s. of Eq. (jL4j) corresponds to 
the free (Hamiltonian) time evolution of the spin, while 
the integral on the r.h.s. describes a dissipative inter- 
action which is nonlocal in time and contains the spin- 
phonon coupling terms. 

The time dependence of the kernel £(t, t 1 ) = — t') is 
evaluated explicitly in Apps. |X| and |b] up to order 0(g 2 ), 
compare Fig. 0b. It depends only on the time differ- 
ence since the Hamiltoni an is time-translationally invari- 
ant. According to Eq. (|A8|), we find that E(i — i') is 
a fast decaying function of time, and we make the sim- 
plifying Markov assumption that p(t) remains essentially 
constant over the time period At within which S(t — i^l 
decays to zero and take p(t) in front of the integral.E3 
This is justified at least for the longest relaxation time 
Ti (see below). For convenience, we also take the upper 
integration limit to infinity. 

After taking p(t) out of the integral in Eq. ( |l4| ) the 
integration over time (up to infinity) can be performed 
and we are left with a constant S. This is evaluated 
in App. [b] and corresponds to the diagrams in Figs. 
and |2l| a and b, see also below. The kinetic equation (|1J) 
becomes 



p(t) = -iL p{t) + S p(t) = Wp(t) 



(15) 



with Lq = [Hs, ■]■ This is similar .to. .the-rnaster equa- 

e^BErta and may be 

(16) 
(17) 



This is simil 
tion formulations in the literatur 
solved for the eigensolutions of W 



i>M{t) = Wp^{t) = \ i pM{t) 
pW(t)=pW(0)-e A **. 



The real parts of the eigenvalues correspond to relaxation 
rates Re{Xi) = — t 4 . 

The eigensolution of W with Ao = corresponds to 
the stationary state, defined as p^ (t) = with p^ (t) = 
p(°\ The diagonal elements of p are the thermal prob- 
abilities to find the spin in the respective states. In the 
m-basis, the static magnetization and susceptibility are 
readily expressed in terms of the diagonal components 



(o) (o) • i i • 
Pin = p m ,m yielding 



M = gpB^2mp$ 



Xo 



dM 



dp 



(0) 



dH z 



(18) 



(19) 



Since the stationary probabilities are given by Boltzmann 
factors we obtain in the absence of tunneling 



Xo 



knT 



^m(m'-m)pV>p<$. (20) 



In the precence of tunneling, this form remains essen- 
tially the same up to transverse fields of the order of IT 



(in this range most of the tunneling splittings are too 
small compared to the level spacing in order to change 
the result considerably). 

For all the other solutions Re(Xi) < and the vectors 
p^\t) correspond to deviations from p(°\ The longest 
relaxation time t± is several orders of magnitude larger 
than T2 and the respective solution is interpreted to cor- 
respond to the over-barrier relaxation. 

Above we have considered the full (reduced) density 
matrix p(t). In appendix ^|, we discuss the choice of the 
basis and argue that the d-basis has certain advantages 
and, e.g., in the case of strong tunneling compard to the 
spin-phonon rates, it allows the restriction to the diago- 
nal elements of p(t)d,d' only. The nondiagonal elements 
become important when the rates for the tunneling and 
the spin-phonon coupling are of the same order of mag- 
nitude. 



C. Ac-susceptibility 

In this section we derive an expression for the dynamic 
susceptibility 



Jo 

dM(t) d(S z (t)) 



(21) 
(22) 



starting from the Hamiltonian Ti and formulating the 
expectation value in Eq. ( |22] ) in terms of diagrams. For 
convenience, we assume in the following that H± can be 
tuned to any (static) value independent of H z and that 
the actual measurement is done by applying a tiny ac- 
excitation field h z (w) on top of the static H z . The more 
general calculation of d M a (t)/d Hp(t'), a,/3 = x,y,z, 
can be carried out along the same lines. 

The derivation with respect to H z (t') acts on the terms 
exp[Ti fdt'H z (t')] oc exp[T*'(-9 Mb) Jdt'S z (t')H z (t'))] in 
Eq. (|) and may occur on either the forward or backward 
propagator of the diagrams (minus and plus signs, respec- 
tively). The derivation takes down a factor Hgp^Szit') 
and we obtain 



X {t-t')=i{gp B Y([S z {t),S z (t')\) 



(23) 



which, when inserted to Eq. (|2l|), is just the Kubo for- 
mula for the linear response to an external magnetic field. 

Susceptibility can be expressed diagrammatically 
starting from either of the equations (^2^ or (^3|). In 
either case, S z (£) is written at the latest point to the 
right. For Eq. (p2[), d /dH z (t') may act on either of the 
propagators, while, for Eq. (p3|), the terms of the commu- 
tator are ordered on the contour reading them from right 
(earlier) to left (later times). It is the most convenient 
to work in the eigenbasis where it is sufficient to consider 
the diagonal elements only, see App. ^|. In this case, if 
S z (t') falls in between S's in the diagram in Fig. Qa, the 
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commutator equals zero - we only need to be concerned 
with the X parts. 

Let us set to = and define E(f",0) as in Fig. 0b but 
with S z (t') (0 < t' < t") inserted onto one of the propa- 
gators. Furthermore, let us assign the part e tuJ (*"-*') f 
the Laplace transform to the definition of Ef t", 0). This 
is then the part of the integrand in Eq. (foVi) appearing 
between times t — and t = t" . On the forward propaga- 
tor, this yields a factor igpB?nd,d' = igpB{d'\S z \d) when 
acting on the state d on the propagator and changing this 
to the state d! (d-basis is not the eigenbasis of S z ). On the 
backward propagator, there is a further minus sign. As 
to the whole diagram, see Fig. [|, the system is in the sta- 
tionary state before to — 0, described by - this cor- 
responds to the requirement that the expectation value 
(S z ) is taken with respect to the stationary state. [One 
could extend the definition of the susceptibility to ac- 
count for experiments with an additional time-dependent 
(e.g. constantly sweeped) magnetic field H z {t) and con- 
sider initial states before t' that are combinations of the 
stationary and the relaxing modes.] After £(t", 0), there 
may be any number of E's between this and the final time 
t - the sum of such series is denoted H(t,t"). Finally 
the resulting diagram/expression is Laplace transformed 
- e l "( t-t ' being divided into two parts, one denoted by 
the dashed line running through the diagram from time 
t" to f and the other one extending from t' to t" and 
already belonging to the above definition of £. 
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Diagrammatic representation of x(w), see text for 



In evaluating the diagram in Fig. ^, we vary t' within 
£(i", 0), and extend the integrations of t and t" to infin- 
ity. This yields 



dte^- t 'H[S z (t),S z (t')])= ]T m d 

d,d',d" 



(24) 



t" 



dt'U(t,t") d ^(t",0) d ^ d ,p^ 



(0) Au(t-t") 

t- 1 



dt" I dt 
10 Jt" JO 

where, for brevity, we denote the diagonal states by just 
one index d, e.g., rrid — rrid.d, see above. The integrations 
can be carried out one by one in the order t' , t (which 
yields U(z = w)), and t" to yield 

= ~(gp B ) 2 (25) 

1 fv ^ (°) 



E 

d,d' .d" 



~iu) + iLo — S(w) 



d,d" 



where S(w) is defined as the Laplace transform of E(f , 0). 
The explicit expressions for H(ui) d d " turn out lengthy 
and are omitted here. 



The resolvent in Eq. ( |25| ) is treated in terms of the 
eigensolutions of W p(t) and the appropriate projections 
of the terms Y,(u>)p^ = ^ i CipW are used. In particular, 
for the component along the relaxing mode p^'(t), the 
resolvent in Eq. (^5|) reads 

1 1 1 



iLo — E(w) — iu> — W(uS) 



-IU) - 



The approximate sign is due to the Markov approxima- 
tion where S(w) ~ S(0), cf. App. ^|. Apart from one n 
this is the well-known factor in expressions for suscepti- 
bility, see e.g. Refs. [39] or [l0| The rest of Eq. ( 25 ) reduces 
to a prefactor roughly proportional to t^ 1 but with weak 
H z and ^-dependences. For the low frequency limit, we 
recover the static susceptibility, i.e., lim w ^ox( CJ ) = Xo 
where \ is given by Eq. (|2ti|). 



IV. RESULTS 

In this section, we concentrate on results obtained 
in the eigenbasis of Hs- This choice is advocated 
by three points. First, it allows us to consider much 
stronger transverse fields (see the discussion on energy 
scales in Part II) than the common approach to calcu- 
late tunnel, splittings in the leading-order perturbation 
theory;oliJ second, the origin of the Lorentzian shape of 
the resonances-, is, seen to arise naturally from the tun- 
nel splittings .EjO The third point is that, in this basis, 
all the relevant properties of the system are captured by 
the diagonal elements p{t)d,d only.L3 This feature also im- 
proves the performance of the numerics. We also present 
results obtained with the full pit) and discuss the differ- 
ences as examples of phonon-induced decoherence. The 
energies are expressed in kelvin, magnetic fields in teslas, 
and rates in s _1 . 



A. Relaxation rates 

The relaxation rate as a function of the longitudinal 
magnetic field H z is shown in Fig. ^. The overall be- 
haviour is of the Arrhenius form, i.e., r ps tq exp\j3U(H z )] 
where To only gives a constant offset to the semilogarith- 
mic figure and U is the effective height of the barrier, see 
Fig. [| On top of this exponential field dependence, there 
are series of resonances located at 



h: 



5Mb 



A + B(m 2 + ml ) 



(26) 



corresponding to the values of the external field that 
brings the states m and m! on-resonance; the reso- 
nances form groups close to the fields H z w nH\ where 
n = m + ml and Hi = see Fig. ^. The broadest of 
the resonances resemble those seen in usual experiments, 
while most of the resonances are too narrow to be seen 
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and, as is discussed below, already the phonon coupling 
is sufficient to suppress them. One of the main topics of 
this paper is the possibility to detect some of the satellite 
peaks by the application of a relatively strong transverse 
field H± as well as to suppress the already visible peaks 
by pointing H± along one of the hard axes. 
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FIG. 9. Relaxation rate r = 1/r as a function of the 
longitudinal magnetic field: without tunneling (dashed) and 
with tunneling (solid). Here, only the diagonal density matrix 
elements are used. The large offset between the two curves 
is due to the broad Lorentzian shapes under the even reso- 
nances: there is a strong direct coupling between the states 
m = — m' = —2 around H z = 0, m = —3, m' = 1 under 
the second series of peaks (H z » 2H\), and m — —4, m! = 
under the fourth series of peaks (H z ~ 4Hi). The sharp res- 
onances arise due to tunnel coupling of lower-lying states m; 
their height reflects the energies E m via the Arrhenius law 
r oc exp[— P(Em — E-io]), and the peak width is given by the 
tunnel splitting. Here the temperature is 2.5K. 

The whole curve in Fig. ^ can be understood in terms 
of the expression 



T = J24 n) Dn(H z )e^ H * 



(27) 



where the index n enumerates the resonances (pairs of 
states m n and m' n ). As discussed in App. there is 
a certain Tq for each resonance; the function D n {H z ) 
is a tunneling-induced Lorentzian that yields the peak 
shapes, and U n (H z ) is the effective barrier height for the 
nth resonance, i.e., U n = E„ ln — .ELio (-E-io is the en- 
ergy of the metastable minimum). The widths of the 
Lorentzian peaks are given by 



5H Z = 



4|A„| 



gm\ m r, 



(28) 



where 2|A„| is the tunnel splitting between the resonant 
states, cf. App. || As is seen in the figure and also sug- 
gested by the inset in Fig. |], the background is not given 



by the over-barrier relaxation but by the "leakage" due 
to the direct coupling of states with to — to' = ±4, e.g., 
to = ±2, to' = =f2 close to H z = and m = ±1, to' = =f3 
close to H z = 2H\. This effectively lowers the barrier 
height to U rj E-2 — E-io- The rates that determine 
Tq , depend only weakly on temperature, the resonant 
states, and H z , cf. App. pi This is in line with the exper- 
imental observation that the estimated prefactor of the 
Arrhenius law, To, varies within an order of magnitude 
for different samples and resonances. 

The sharpness of the narrowest peaks in Fig. ^ is an 
artefact of the reduced model used so far, i.e., neglect of 
the nondiagonal matrix elements, and let us next con- 
sider two things affecting these peaks: decoherence due 
to spin-phonon coupling and the broadening effect of 
a transversal magnetic field Hj_. For other contribu- 
tions sucli as the hyperfine and dipolar interactions, see 
Sec.lfvC. 



The decay of the nondiagonal elements of the density 
matrix is a classical definition of decoherence and this is 
also what is seen here when the above calculation is per- 
formed using the full pit). Inclusion of the off-diagonal 
elements still produces Lorentzian peaks but now the 
narrow resonances - with widths of the same order of 
magnitude as the spin-phonon coupling - are broadened 
and reduced in height. The narrowest peaks may even be 
completely suppressed, see the solid-line curves in Fig. 
Despite its intuitive appeal, the suppression of the narrow 
peaks should be taken only qualitatively because of the 
limitations of the Born approximation used in treating 
the spin-phonon coupling. For this reason, the following 
considerations are restricted to regimes where the effect 
of the nondiagonal elements of p(t) is negligible. 
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FIG. 10. Magnification of the first and second clusters of 
peaks of Fig. [] magnified (solid line). The different curves 
correspond to: H± = 0.0T (solid), 0.01T (dotted), 0.05T 
(dashed), 0.1T (long dashed), and 0.2T (dot dashed); = 0° 
for all the curves. The peaks further to the right correspond 
to resonances between lower states, hence the higher maxima. 
The peaks (from left to right) correspond to the resonances 
m — 1, .., 4, m! — —m — 1 in the left figure (the two first ones 
are merged together) and m = 2, .., 5,m' = — m — 2 in the 
right figure. Here we used also the offdiagonal elements of p(t) 
with the result that some of the sharper peaks are reduced 
in height or even suppressed due to spin-phonon coupling. 
The peak widths show strong iix-dependence except for the 
second and fourth peaks on the right - these arise from B\ and 
Bf-type of coupling. Once the peak widths due to the tunnel 
splittings exceed the phonon-induced width, the peak heights 
are essentially determined by the energy U n = E n — E-io 
needed to reach the resonant states. 

The suppressing effect of the spin-phonon coupling be- 
ing essentially a constant, the sharper resonances can be 
made observable by broadening them with a transverse 
field H±; even grouncUstate tunneling has been observed 
in high enough fields J13 The various tunnel splittings are 
shown in Fig. || as a function of H± pointed to three dif- 
ferent directions <j>. It can be seen that for small values 
of the field, the splittings are essentially independent of 
the chosen angle, see also Fig. [Lll Furthermore, tunnel 
splittings between states that can be coupled with sole 
Bi terms of TIt are quite insensitive to the transverse 
field below 0.1-0.2T. Figure |o| illustrates the effect of 
H± = H x onto the two first series of peaks increasing H x 
monotonously increases the tunnel splittings. The reso- 
nances that initially (H x = 0) were broader than Tq 1 ^ 
are seen to retain their height; the narrower peaks, such 
as the two shown in the insets, are strongly broadened 
with increasing H x and also their heights increase once 
their widths become larger than the spin-phonon rates. 

A transverse field cannot only increase the tunnel split- 
tings but it can also reduce them, see Fig. || In the ideal 
case, the splittings can even be suppressed by varying the 
angle <f> where the transverse field is pointed. For the ze- 
roth peak, i.e., for H z ks 0T, all ten resonances occur at 
the same position and the suppression of any one of them 
is obscured in the relaxation rate curves by the simulta- 
neous broadening of other resonances. For this reason, 
let us first focus on the resonances occuring at finite H z 
and return to the H z w case in the next section. 

Figure pd| shows the tunnel splittings A 4j _ 5 and A 3 ,_ 5 , 
respectively, as a function of H± and for four different 
values of <j). The former corresponds to the resonance 
shown in the left inset of Fig. |l^, while the latter corre- 
sponds to the broad resonance in the right panel of Fig. [h] 
which is also the one seen in experiments at H z w 1H\. 
Figure |ll| also exemplifies some general properties of the 
tunnel splittings. As already mentioned above, the A's 
depend only weakly on <f> for low H±, while for larger 
fields and <f) close to the hard axes, the tunneling am- 



plitudes are successively increased and reduced as H± 
is increased. On the other hand, with H± held fixed 
to one of the suppression points, the A's oscillate as a 
function of <f>. Similar to recent work on Fes, see e.g. 
Ref. [li|, these phenomena are attributed to alternating 
constructive and destructive interference.|G>f the geomet- 
rical phase in the tunneling amplitudes.! 45 }! 46 ! The number 
of minima for a given A m , m > appears to be, as was also 
pointed out in Ref. given by the number of i?4-terms 
involved in coupling the states m and to'. This ranges 
from zero at the top of the barrier to five for the gound 
states to = ±10, cf. Fig. f|. 
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FIG. 11. The tunnel splittings A 4 ,_ 5 (left) and A 3 ,_ 5 
(right) as a function of the transverse magnetic field when H z 
is kept fixed at the position of the corresponding resonance, 
0.4581T and 0.9045T, respectively. In both figures, the four 
curves correspond to different angles <j>: 0° - dot-dashed, 40° 
- dotted, 43° - dasheded, and 45° - solid; the vertical axes 
are the same. 



Figure 12 shows the relaxation rates corresponding to 
the resonances of Fig. [ll] with _ffj_'s chosen to match the 
first points of suppression H™< m , The resonance width is 
found to be very sensitive to <f> quite as expected and, for 
<f> = 45°, the left-most peak becomes suppressed. Note, 
that the suppression is also quite sensitive to the value 
of H± as well and the higher peaks in Fig. 12 are hardly 
affected at all by changes in (p. 
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FIG. 12. Relaxation rate r _1 as in Fig. [To] but now for 
stronger transverse field H± and for different angles <f>. The 
line types and angles are as in Fig. [ll] and the values of 
H± = 0.4311T and 0.2643T are chosen to match the first 
suppression points (for the actual suppression occurring at 
the minimum value for A the values for H± should be given 
with even higher accuracy). Temperature is 2.5K. 



B. Ac-susceptibility 

In explaining experimental results, the susceptibility 
x(w) is often described by the formula 



x(H,T,4>,e,u) = 



Xo(H,T,( 



1 — ilut(H, T, i 



(29) 



where t{H, T, (f>, 9) is the relaxation time, u> the frequency 
of the small excitation field, and xo(H,T,9) the static 
susceptibility, cf. Ref. [|. We emphasize here the in- 
dependence as we wish to investigate the interference 
effects and resulting oscillations in t\ (below we do not 
explicitly write the parametric dependences of r but they 
are implicitl y assu med). Equation (^9|) relates to the re- 
sults of Sec. [II C as follows. Equation (p5| ) formally ac- 
counts for all the modes i and for each mode divides into 
two contributions, one being 1/(1 — MJTj) and the other 
the respective prefactor. It turns out that, for most pa- 
rameter values, the over-barrier relaxation, i.e., the mode 
z = l strongly dominates over the others and one can use 
Eq. (H|) to a good approximation. Figure [l3| illustrates 
the two contributions to Eq. (|2^); they form the basis for 
understanding the following results. In all figures, the 
susceptibility is expressed for a single spin and in units 
of KT -2 . 

In the present model, the main correction to Eq. j29| ) 
corresponds to an intra-valley mode describing transi- 
tions between the two lowest states on the same side of 
the anisotropy barrier. This correction increases with 
increasing H z and temperature and, for the parameters 
and the curves of xo m Fig. |l3|, becomes of the same 
order of magnitude as the actual relaxing mode once 



H z ss 0.9— IT. However, such a mode corresponds to very 
high frequencies and is neglected below. There could be 
also other sources of high-frequency contributions such as 
dipole-dipole flip-flops, nuclear spin dynamics, or moving 
impurities but these are beyond the scope of this work. 
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FIG. 13. Static susceptibility \o as a function of the lon- 
gitudinal magnetic field shown for different temperatures: 2K 
(solid), 2.5K (dotted), 3K (dashed), 4K (long dashed), and 5K 
(dot-dashed). The dynamic susceptibility x(lo) is essentially 
composed of \o multiplied by the real and imaginary parts of 
1/(1 — icor) shown in the inset as the solid and dashed lines, 
respectively. 

As to the actual results, it is in principle sufficient to 
combine the curves for r of the previous section with 
those shown in Fig. [l3| The shapes of the real and imag- 
inary parts of 1/(1 — iujt) readily imply how the resulting 
susceptibility behaves if one fixes to - it turns out that 
all the structure found in t _1 in the previous section can 
also be found in x(o;). First, the response is the most 
sensitive to changes in r when lu « r _1 or lot w 1, cf. 
the inset of Fig. [13[ Second, Ke(x(H z ; ui)) replicates the 
shape of t(H z ) and exhibits peaks and valleys as r _1 
increases and decreases. This is the case with the imag- 
inary part, lm(x(H z ; lu)), as well, but only as long as 
lot > 1, i.e., as long as we remain on the right hand side 
of the peak in Im{x{H Z) m)) . For lower u>, lot can cross 
the maximum point and the picture with the peaks and 
valleys turns upside down. These points arc illustrated 
in Fig . |l4| for a case corresponding to one of the curves 
in V 
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FIG. 14. The real and imaginary parts of susceptibility 
for the same parameters as the dotted curves in the left pan- 
els of Fig. i.e., H± = 0.4311T, cf> = 40°, and k B T = 2.5K. 
The curves correspond to frequencies = 0.02Hz - \' (dot- 
ted), x" (dashed) - and LU2 = 0.5Hz - \' (dot-dashed), %" 
(solid). The frequencies were chosen such that wi < r _1 for 
the both of the peaks in Fig. ^ - hence the notches in x", 
while u>2 > i"" 1 for all H z . 

Let us next consider a different scheme, keeping H z 
and 4> fixed and varying H± starting from zero field. The 
case H z w 0T is of particular interest because it allows 
comparison to recent experiments by Wernsdorfer done 
on the related material Fes, cf. Ref. [|7|- the experimen- 
tal data showed clear oscillations of the relaxation rates 
as a function of H±_ in the thermally activated regime. 
By choosing parameters in the feasible range of these 
experiments, we find somewhat similar oscillations also 



for Mni2. Figure 15 shows first the relaxation rates for 
two different combinations of temperature and the lon- 
gitudinal field: T = 3. OK and H z = 0.2mT to show the 
behaviour close to (or at) the peak maximum and at a 
lower temperature, and T = 5. OK and H z = lOmT as 
an example of the behaviour further away from the max- 
imum and at a higher temperature. Also the <f> depen- 
dence is shown. We would like to stress the novelty of 
this result: such oscillations have neither been observed 
nor predicted before. 
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FIG. 15. Oscillations in the relaxation rate as a function 
of H± . The lower and upper groups of curves correspond to 
H z = 0.0002T, k B T = 3.0K and H z = 0.01T, k B T = 5.0K, 
respectively. The individual curves correspond to different 
angles 0: 0° (dashed and dot-dashed), 43° (dotted), and 45° 
(solid and long dashed). 

For <fi — 0°, the gentle oscillations resembling a 
strongly-smeared staircase stem from the fact that the 
lower (in terms of energy; higher in terms of rates) reso- 
nances are broadened and one by one start to dominate 
the relaxation, see the dot-dashed curves in Fig. [tf] for 
the tunnel splittings and Fig. [l?] for the general idea. All 
the additional structure, seen for <fi — 45°, is due to the 
oscillations in the tunnel splittings and the suppression 
of some of the resonances. The positions of the notches 
can be compared with the structure of A TO) _ m in Fig. ^ 
and one finds that for smaller H z the relaxation takes 
place via the lower-lying resonances such as m = ±6 and 
m = ±7; further away from the maximum, the broader 
resonances between states m = ±4 and m = ±5 act as 
the dominant relaxation paths. 
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FIG. 16. Tunnel splittings A m> _ m for four resonances 
relevant for Fig. |l5|. The dot-dashed and solid lines denote 
the angles (f> — 0° and 45°, respectively. 
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FIG. 17. Log-log scale schematic of Lorentzian curves on 
top of each other similar to the H z ~ 0T situation for Mni2. 
The higher the curves are the narrower they get and, on the 
other hand, the higher the H z the lower the observed reso- 
nances are. The relaxation rate is determined by the fastest 
possible rate, depicted in the figure with the thick solid line. 



Figures |l8| and 19 show the real part of the suscepti- 
bility x' corresponding to the two cases of Fig. |l5[ The 
purpose of the different frequencies is to show that also 
here one can choose the structure of interest and study it 
by tuning the frequency to fulfill lot k, 1. In both figures, 
there are regimes where the susceptibility can be varied 
by a factor of five; in Ref. [l7| the oscillations were quite 
clear already with the amplitude being a mere 20% of the 
signal. 
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FIG. 18. Real part of the susceptibility for H z = 0.2mT, 
k B T = 3. OK, and for the angles <j> = 0° (dashed), 43° (dotted), 
and 45° (solid). The upper three curves correspond to the 
frequency lu = 10Hz and the lower ones to lu = 50Hz. The 
line types correspond to those for r _1 in Fig. 
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FIG. 19. Real part of the susceptibility for H z — lOmT, 
feT = 5. OK, and for the angles (f> = 0° (dot-dashed) and 
45° (long dashed), i.e., the line types correspond to those in 
Fig. ttq. The upper pair of curves corresponds to lu = 250Hz 
and the lower one to lu = 1250Hz. 



C. Discussion — relevance to experiments 

So far we have considered the simple model compris- 
ing a single spin coupled to a phonon bath but, as was 
pointed out in the introduction, in real samples there 
are also other kinds of interactions. In this section, wc 
consider the additional features arising from the hyper- 
fine and/or dipolar interactions and aim to point out the 
experimentally relevant aspects of the results obtained 
above. 

Let us first recall some experimental facts concern- 
ing relaxation measurements and results - these under- 
lie also the understanding and appreciation of the sus- 
ceptibility measurements. Relaxation rates are typically 
measured by first magnetizing the sample to satura- 
tion, and then reversing the direction of the field and 
measuring the resulting magnetization as a function of 
time. The initial relaxation is observed to be nonexpo- 
nential - this is attributed to dipolar interactions, sea 
below - while, at later times, it becomes exponential.E£l 
Several authors have proposed an extended exponential 
M(t) = M(0) exp^t/r)' 3 ] to account for both of these 
regimes with just one additional fitting parameter [3. In 
Ref. [is] it was found that (3 varies from /3 ~ 0.5 below 
2. OK to (3 k 1 - usual exponential relaxation - roughly 
above 2.4K. The thus obtained relaxation rates exhibit 
a series of broad Lorentzian-shaped resonances; their 



12 



height and location correspond to tunneling-assisted re- 
laxation 3-4 levels below the top of the barrier. 

This shows two clear differences as compared to the 
present work: in experiment, the relaxation may be non- 
exponential even though the single-spin model always 
yields exponential behaviour, and no satellite peaks are 
observed (see Ref. ^5] for exceptions) . In order to under- 
stand these discepancies, let us first consider the effect 
of the nuclei via the hyperfine interactions and then the 
intermolecular dipolar interactions. 

Hyperfine interactions. In Mn 12 , all the manganese 
nuclei have magnetic momenta and the hyperfine inter- 
action between the nuclei and the molecular spin state is 
relatively large, of the order of lOmT. Recently, several 
authors have investigat cd. . h tj w . rtteh affects tunneling and 
the relaxation in Mni2.tHl"ii3'cJS For the present pur- 
poses, the relevant effect of the hyperfine interactions is 
to induce an intrinsic Gaussian broadening, of the width 
CThyp w 6mT,tJ to all levels, i.e., the nuclear spins are im- 
portantly dynamic and their influence on the molecular 
spin cannot be reduced to a rigid but spatially varying 
background field. Simultaneously with the broadening, 
the hyperfine interactions reduce the tunneling ampli- 
tudes of the resonances for which A < <7hyp ~~ this should 
lead to reduced peak heights in t~ x (H z ). 

With this in mind, let us consider the different regimes 
in terms of the relative magnitudes of the tunnel splitting 
and the hyperfine broadening. First, for resonances with 
A 3> 0hyp the shapes of the resonances are expected to 
be Lorentzian with the widths determined by the tun- 
nel splittings A. This is the regime, where all our re- 
sults apply. In the other extreme, A <C Chyp, the res- 
onances should be essentially suppressed providing one 
possible explanation why the sharp satellite peaks are 
not observed in experiments. Note that the minuscule 
phonon-induced broadening or dephasing is hidden un- 
der the hyperfine broadening and cannot be seen. In this 
regime, one could in principle try and extend the present 
theory by adding by hand a strong dephasing term to 
the nondiagonal density matrix elements. The interme- 
diate regime where A w <Xhyp is the most interesting of 
the three cases. In this regime, the peak shape should be 
a combination of Lorentzian and Gaussian curves and, 
depending on which one of A or <7hyp is larger, one of 
the shapes should dominate. In the tail region, i.e., away 
from the peak maxima, the Lorentzian tails dominate 
and it has been suggested that this together with exper- 
imental error bars may obscure the resolution between 
the two types of curves, cf. e.g. Ref. ^3|. 

The immediate conclusion from these considerations 
is that, if the application of H± broadens some of the 
tunnel splittings A m , m i to exceed Oh yp , the correspond- 
ing resonance should become observable. If, on the other 
hand, the transverse field is applied along one of the hard 
axes, a given resonance becomes suppressed for certain 
special values of H±; in the presence of the hyperfine 
interactions this should happen already when A m , m < be- 
comes smaller than Ohyp. The intermediate regime can 



be intentionally achieved by tuning the tunnel splitting 
from being well below <7h yp to above it. This may pro- 
vide means to probe the Gaussian broadening, see also 
the subsection below focusing on the advantages of sus- 
ceptibility measurements. 

Dipolar interactions. The intermolecular spin-spin 
interactions are of dipolar form and they are weaker in 
Mni2 than, e.g., in Fes- Due to their short range, the 
dipolar fields can vary in space changing the local field 
at the position of the individual molecules. In our view, 
the essential difference between the hyperfine and dipolar 
interactions can be stated as follows: even if one could 
measure the response of a single molecule, this would al- 
ways be dressed by the level broadening and reduction 
in tunneling amplitudes due to hyperfine interactions in- 
trinsic to each molecule; the dipolar fields, on the other 
hand, just change the molecule's local electromagnetic 
environment. 

In experiment, the relaxation of the magnetization 
M(t) leads to time-dependent dipolar fields and, in order 
to describe the relaxation correctly, it would be necessary 
to solve for M(t) self-consistantly, for simulations see e.g. 
Refs. [19] and |2(]. However, it is this time-dependent field 
that provides an explanation to the initial nonexponen- 
tial relaxation. For example in Ref. |l5|, it is therefore 
concluded that the deviation from a single-exponential 
relaxation, i.e., f3 =/= 1, demonstrates the important role 
of the dipolar interactions and dynamics of the spin dis- 
tribution. It should be kept in mind, though, that also a 
static distribution of local fields (be it dipolar or not) - 
and hence relaxation rates r _1 (i?' ocal ) - leads to a super- 
position of exponential rates, which looks nonexponen- 
tial. Such a distribution of fields also hides all features 
in t (H z ) that are sharper than this distribution. 

The time-dependence of the dipolar fields owes to the 
fact that the sample is first magnetized and, as the field 
direction is abruptly reversed, the dipolar distribution 
finds itself far from equilibrium and quickly starts to re- 
lax. The reason for such experiments is the strong re- 
sponse from almost all the spins. 

In anticipation of the discussion on susceptibility, let us 
consider the dipolar distributions at equilibrium. By dis- 
tribution we mean spatial variations in the dipolar field at 
the locations of the individual molecules. First, for H z w 
OT, the annealed (not quenched) distribution is random 
but due to the low temperatures, k^T <C E±g — E±io, 
almost all the spins are aligned with the easy z-axis - ran- 
domly pointing to the positive and negative directions - 
and only contribute to the local longitudinal field. On 
the other hand, the equilibrium magnetization is close to 
saturation already for H z « H± , thus drastically narrow- 
ing the dipolar distribution - e.g. for T = 3. OK (5. OK) 
and H z = Hi, more than 95% (85%) of all the spins 
are aligned parallel to H z and all the molecules feel es- 
sentially the same field. Such distributions have been 
experimentally verified in Fe§ , cf. Ref. 

Susceptibility. The influence of the dipolar dynam- 
ics on relaxation can be avoided almost completely by 
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measuring linear response to a small ac- field, i.e., the ac- 
susceptibility, instead of M(t). This has the advantage 
that the system is probed in its equilibrium state and 
ideally by a small enough field that in itself does not per- 
turb the equilibrium. Therefore we propose that the sus- 
ceptibility measurements provide a gentle or noninvasive 
means to probe the relaxation dynamics in absence of the 
time-dependent dipolar distribution. Furthermore, while 
the hyperfine interactions cannot be tuned, the static dis- 
tribution can be made markedly narrower by a finite H z , 
see above. As on the other hand xo decreases with in- 
creasing H z , the first group of resonances, close to Hi, is 
especially attractive for investigating the oscillations in 
the relaxation rates as well as the hyperfine fields them- 
selves. All of the resonances around Hi depend strongly 
on H± and can be broadened such that A > Ohyp making 
them observable; the peaks can also be selectively sup- 
pressed if <f> ~ 45°. The hyperfine fields may even sim- 
plify the observation of the suppressions as very narrow 
peaks are strongly reduced in height. For a sharp dipolar 
distribution, we expect that also the crossover between 
Lorentzian and Gaussian shapes of r _1 (iJ z ) should be 
observable when changing the tunnel splittings with the 
transverse field. 



tain values of H±. This phenomenon attributed to the 
interference of the geometrical or Berry phase of alter- 
native tunneling paths, with a destructive interference 
leading to the suppressions. As the second major result, 
we predict that these oscillations in the tunnel splittings 
should be observable both in the relaxation rates and the 
susceptibility. 
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V. CONCLUSIONS 

To conclude, we present a diagrammatic description of 
the spin dynamics of the molecular magnet Mni2- The 
work focuses on the regime of thermally-activated tunnel- 
ing, i.e., T > 2. OK, and emphasizes the phenomena that 
could be observed for strong transverse magnetic fields. 
In the calculations, we study the dynamics of a single 
spin S = 10 coupled to a phonon bath. The role of the 
phonons is in the thermal activation of the spins to states 
with higher energies and larger tunneling amplitudes. 

As the first main result, we calculate the dynamic 
susceptibility starting from the same microscopic 

Hamiltonian as is used for the relaxation rates. Sus- 
ceptibility is found to reflect the rich structure found in 
r -1 (-ff z ) and we argue that susceptibility measurements 
are in fact more sensitive and better controlled in terms 
of time scales and also the dipolar interactions than the 
relaxation experiments. 

All the results obtained are calculated using the eigen- 
basis of the spin Hamiltonian, which naturally accounts 
for strong transverse magnetic fields. A strong transverse 
magnetic field enhances tunneling through the anisotropy 
barrier and enables relaxation via eigenstates further 
away from the top of the barrier. In relaxation or sus- 
ceptibility measurements, this would lead to shifted and 
higher resonances. The tunnel splittings are found to 
be very sensitive to the azimuth angle <fi of the trans- 
verse field H±. It is found that, in the directions <f> — 
7r(2n+l)/4, the tunnel splittings exhibit alternating min- 
ima and maxima and become totally suppressed at cer- 



APPENDIX A: DIAGRAMMATIC RULES 

In this Appendix, we list the rules for evaluating the 
diagrammatic expressions for the kernels £ = J dtT,(t) 
arising in the evaluation of (S z (t)) and the kinetic equa- 
tion. We assume that the Markov approximation is valid 
and check the results for self-consistency. 



1. Real-time representation 

In the time domain, the rules for £(t) are 

1. Propagators 

Assign a factor exp[— i7is{U+i —ti)] to each piece of 
a forward propagator between two vertices at times 
U and tj+i (ti < ti + \). For a particular diagram 
with specific states at the ends of the propagators, 
cf. Fig. pOja, take the element 



iH s (U 



(Al) 



For a backward propagator, change the sign in the 
exponent and invert the order of the states uij and 
m\. In the m-basis, the tunneling contained in TLs 
changes the states along the propagator but in the 
d-basis the above exponential reads 



iE d .(ti 



-U) 



(A2) 



A compact way to account for the tunneling in 
Eq. (Al) is to rewrite it with the help of Eq. ( A2) 



as 
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-iE d At i+ i-tt) 



(di\m,i). 



(A3) 



In so doing, wc come to incorporate the tunneling 
elements exactly to the corresponding diagram. 

Vertices 

Assign a prefactor —i (+i) to each vertex lying on 
the forward (backward) propagator. For each pair 
of vertices at times and t i+1 and coupled by a 
phonon line, cf. Fig. pO|b, assign a factor 



G 



mm ,m\m\ 



\J $m , m ' $m i , rn i * ^-^,^ / j 



(A4) 



where £ 



and £ = mi — mi, see App. |§| 



3. Spin-phonon interaction lines 

A phonon line connecting a pair of vertices corre- 
sponds to a phonon correlation function - this is 
worked out in detail in App. [b]. It is more conve- 
nient to first calculate this in the energy represen- 
tation 



A 2 



12pc 5 ?i 4 



a f3ui 



and then Fourier transform it 

poo 

T(t-t')= / duj T(uj) e M*-*') 



A 2 



12pc 5 h 4 - 1 



(A5) 

(A6) 
■^ t ' t '\ (A7) 

Here A is the same as in TL Z , cf. Ref. |l8|, p = 
1.83 • 10 3 kg/nr is the density of Mni2, c is the 
sound velocity (the only "fitting" parameter of the 
theory), and u> is the phonon energy, lo > mean- 
ing absorption and uo < emission. Note that the 
prefactor here differs from that found in the litera- 
ture for the spin-phonon rates. This is just because 
for convenience we assign part of it to the vertices, 
cf. Eqs. (B8)-(B7). Fourier transforming the cor- 



relator into the time domain gives 



T(at) 



A 2 



Ylp&Ti 



(A8) 



- 1 ' i 



1 + 2 ch- 



( T<Jt \ 
\ P J 



sh' 



( E£l\ 

\ a J 



7T<5 W ((rt)} 



where a = +1 (-1) if the vertex with the earlier 
time lies on the forward (backward) propagator; 
the function 5"'(t) is the third time derivative of the 
delta function. This expression diverges at t = 
but this is just an artefact that is removed by in- 
cluding a high-energy cutoff in Eq. (|A7|). For in- 
stance, an exponential cutoff e^l D would lead to 



the replacement t — > t — if D in Eq. (A8) and thus 
avoiding the divergence. We are not interested in 
the limit t — > but rather the longer-time be- 
haviour of the correlations: the correlator decays 
exponentially on the time scale of l/k^T which is 
5-10 orders of magnitude faster than typical relax- 
ation times r of the reduced density matrix. This 
observation is used to justify the Markov approxi- 
mation in the text. 

4. Summations and integrations 

Sum over the states internal to the diagram. For 
example, in evaluating the element S mmi)7n / m ^ in 
Fig. |2fj|b, sum over the states m and mi, and in- 
tegrate over all the times (time differences) in the 
diagram. 

We might add a fifth rule to capture a general prop- 
erty of combinations of diagrams: a series of irreducible 
diagrams is evaluated iteratively such that all the infor- 
mation needed from earlier processes/diagrams is incor- 
porated into the reduced density matrix p(t). As an ex- 
ample, the diagram in Fig. p0| b equals 

m,mi 

• *^^(m'i\di) (di\mi) (m\d) (d\fh) 

dt , e -i(E d -E dl )(t-t') r ^ _ t ) p (t') m , mli . 



a) 



b) 



> • • 



ti+l 



mm m 
> i , — 



— «g- 

m] 



mj mi 



FIG. 20. Sample diagram for illustrating the rules above. 
The vertical dashed lines in (b) are guides for the eye: the 
rules in the text are used to evaluate the piece of the diagram 
between times t' and t that contributes to E(t' — t). 



2. Energy representation 

Once the rules for calculating £ have been laid down 
in time domain, it is more convenient to shift into the 
energy representation with respect to the phonon ener- 
gies and to do this in the <i-basis. The rule 2. remains 
as it is, Eq. (A5) in the rule 3. is already in the energy 
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representation, and we only need to reformulate the first 
and fourth rules. 

1. Propagators 

Assign a resolvent 

(A10) 



(juj - E d i + E d 



to each piece of a diagram temporally between two 
vertices, i.e., irrespective of the propagator they lie 
on. Here to is the phonon energy, and Ed' and Ed 
are the energies of the states on the forward and 
backward propagators, respectively; a = +1 (—1) if 
the vertex with the earlier time lies on the forward 
(backward) propagator. The factor r\ > arises 
from the adiabatic turning on of the interaction and 
is taken to zero at the end. 



a) 
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d'i 
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d i 



d j 



d j 



FIG. 21. Diagrams may sometimes be pairwise combined 
to simplify the calculations. 



4. Summations and integrations 

Sum over all internal indices and integrate over u>. 
Due to the factor irj, the integrations are to be 
understood as combinations of Cauchy's principal 
value integrals and delta functions. 



In diagrams such as those in Fig. |21| a the prefactors due 
to the vertices are equal. Therefore these can be pairwise 
combined and their integral parts written together as 



dus T(oj) 



(All) 



-uj + (E d ' -E dl ) + ir) u)- (E d ' x - E d ) + irj 



If d' = di and d = d\, i.e., if the density matrix is di- 
agonal before and after the transition, this expression 
simplifies into 



7T / duT(w)S(u-(E d >-E d )) (A12) 
irA 2 AE 3 



nT(E d ,-E d ) = 



12pc 5 h 4 ' eP AE -r 



In the last step we inserted the definition of T(uj) from 
the equation (B6) and identified AE = E d i — E d as the 
energy required for th e tra nsition from the state d to d' . 
This expression, Eq. (A12), together with the prefactor 
Gd'd^id^ is the phonon induced transition rate Y, d ,d' used 
in the litera ture and in Par ts HI and IV of this paper. 
From Eqs. ( All ) and ( Ml ) we see that the energy is 
only conserved (the delta function, cf. Fermi's golden 
rule) under some special circumstances and, in general, 
the rate need not be the simple real function used in the 
literature. 



APPENDIX B: THE SPIN-PHONON RATES 

In this Appendix, we outline the calculation of the 
phonon-induced transition rates £ between different spin 
states. In order to calculate S in lowest order in the 
spin- phonon coupling constants gt, we need to evaluate 
the contractions (H S p(t)H S p(t'))ph, where the expecta- 
tion value is taken with respect to the phonon degrees 
of freedom. Let us first insert Eq. (g) into Eq. (f|) and 
explicitly write down all the resulting terms 



n sp (t) = iJ2 



[bl e 



ka 



l itor t~\ ik-r 
her H e 



(Bl) 



■{gx[Sl(t)-Sl(t)][e^k x -e^k y ] 
+ il{ Sx (t),S y (t)}[e^k y + e^k x } 

+ 9 -l({S x (t),S z (i)}[^ ] k z +e^k x ] 

+{S y (t),S z (t)}[e^k z + e<fh v }) 
+ °±({S x (t),S z (t)}[e { ° ) k z -e < ?h x ] 

+{S y (t),S z (t)}[e y °h z -e z °h y ])}. 

In evaluating the contraction, we need to consider all the 
possible states before and after the action of Tt sp (t^) as 
well as all the possible orderings of the times t and t' 
along the contour. As a characteristic example, let us 
consider the diagram in Fig. |2p| b and hold to the same 
usage of indices. The time dependences of the spin oper- 
ators are accounted for by the propagators and they can 
be put aside for the moment. The remaining part with 
the spin operators is fully determined by the spin states 
before and after the vertices yielding 



(mi 



m - m' ^mi ,mi • 



(B2) 

Here £ = m — m', £' = mi — fh\ can take values ±1 or ±2 
due to the structure of H sp , and, for brevity, S ,± '^' = S± . 
The quantities on the right hand side of the equation are 



1G 



Si) 



(2m' + y/S(S+l)-m'(m'+£) (B3) 



for £ = ±1 or 



Srn, m - = {[S(S+1) -m'[m' + v)] 

[S(S + 1) - {m! + v){m' + 2v)}} 



(B4) 



1/2 



for £ = ±2, v = sign(£). 

The part of the contraction depending on the phonon 
degrees of freedom is evaluated assuming acoustic 
phonons with a linear dispersion relation and three modes 
indexed by a = 1,2,3 (one longitudinal and two trans- 
verse modes). The polarization vectors are assumed 
unit vectors. The resulting expression consists of two 
parts. The first one of them contains all the information 
concerning the phonon spectrum, energies, and temper- 
ature, 



/oo 
duj r(w) 
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j(t-t') 



du> 
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12pc 5 n 4 - 1 



Ju(t-t') 



(B5) 



(B6) 



This part is defined such that it is independent of the spin 
states and only contains the term A 2 from the coupling 
constants that turns out to be constant for all the rates. 
In the diagrams, this corresponds to the spin-phonon in- 
teraction line. The second part, 



o, foriei ± in 

(\- = \ 1, £ = £' = I' 



(B7) 



however, does depend on the spin states and implies ad- 
ditional selection rules. It should be noted that the third 
clause allows for £ ^ £' if |£| = |£'| = 2. For convenience, 
we combine the term CV f/ with the spin operators and 
obtain the term 



G 



mm ,mimi 



(B8) 



which in the diagrammatic language corresponds to a 
pair of vertices. 

In the (i-basis, the time dependence of the s pin oper- 
ators is of a simple exponential form, cf. Eq. ( |A3[ ), and 
the diagram in Fig. EOJb can be evaluated to yield 



E(i-t') 



(B9) 



mmi ,m'm x 

E G mm',mimi ' ^{m'^) {d^m^ {m\d) (d\m) 
fh.fhi d,d\ 

e - l(Ed -E d ,)(t- t >) r (t'_ f ). 

The actual transition rates are obtained by integrating 
Eq. (BE) over the time difference r = t — t'. This takes 
us to the energy representation, combining the exponen- 
tial factor in Eq. (fB9| ) with the e 4 "'* _t ' of the Fourier 
transform in Eq. QB6|). The intergration over r yields 



mmi ,771m. 



(BIO) 



r( w ) 



d.di 



■ i I duj - 



-uj - Ed + Ed' + ir] 



As the final s tep i n calculating S, let us evaluate the 
integral in Eq. (BIO) or 



I (T (AE,P)= / dio 



= P I duj- 



(Bll) 



AE 3 



•P" - 1 a(uj - AE) e 0AE _ 1 



for the general case. Here a = ±1. The imaginary 
part of the integral (the last term) is, up to a pref- 
actor, just — z7rr(A_E) and is independent of a. The 
real part in turn can be evaluated using the calculus 
of residues. By combining the integral along the real 
axis with an infinite semicircle in the upper half plane 
to form a closed contour we obtain Ke{I a (AE, (3)} — 
— 27rlm ^{enclosed residues}. 

The real part of the integral I a (AE,(3) is divergent 
due to the w 3 -term and some kind of a cutoff procedure 
is needed here. We have chosen a functional cutoff using 
instead of I(AE, 0) the integral 



I a {AE, f3,D) 



duj- 



D 



(B12) 

2 x 2 



a(uj - AE) + ir) \D 



i.e., the cutoff function is a Lorentzian squared and the 
additional parameter D is the cuto ff parameter of the 
Lorentzian. The integrand in (B12) has single poles at 
lu = AE — iam and at to = 2m7ri//3, where m is a positive 
integer; there are also second-order poles at lj = ±iD. 
The residues from the first pole are real and can be 
neglected. Collecting the other poles in the upper half 
plane and evaluating the residues yields 



Re {^(AS, /?,£>)} = -2tt 
where 



m—l 



cB = 



Cr(3) 



J3 2 AED 4 

(2tt) 3 ' ( AE(3 
\ 2tt 
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(B14) 



4(A£ 2 + D 2 ) 



(3AE 



4 sin 2 (B°) 
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If D is accidentally chosen to equal 2irm* / '/3 for some m*, 
the residues C^l + C' 3 ' should be replaced by 



C* =a 



D 3 3AE 5 - UAE 3 D 2 - AED 4 
16/3 



(AE 2 + D 2 ) 3 



(B16) 



In the d-basis and using only the diagonal elements of 
the reduced density matrix, the integrals I a (AE, (3, D) 
can always be combined such that their real parts can- 
cel each other. When using the nondiagonal elements of 
p(t)d,d', as well, this is no longer true and the real parts 
of I a (AE, fj, D) give rise to D-dependent shifts in the 
energies Ed- In this work, we have used D's of the order 
of the anisotropy barrier, i.e., 50-100K, and found that 
this yields tiny shifts also in the relaxation rate curves 
but leaves the qualitative picture unchanged. In the fig- 
ures in the section IV, I? is chosen as low as 25K in order 



to simplify the comparison between the calculations with 
and without the nondiagonal elements. 

With all the contributions to £ written down, we can 
find an estimate for the order of magnitude of the ele- 
ments (the individual rates) S mmi t m'm' ■ The most in- 
teresting piece of information for each state is the largest 
rate coupling that state to other states - this rate plays 
a key role in justifying the neglect of the nondiagonal 
states in p(t), see App. ^|, as well as in the suppression 
of the narrow resonances found in the text. 

The prefactor in T(ui), cf. Eq. (B6), amounts to 
7.0 • 10 5 c _5 s 5 m _5 K" 2 , where the sound velocity c is ex- 
pressed in meters per second. The units are chosen such 
that, when lu is expressed in kelvin, also T(uj) is given 
in kelvin. For u> > (and also u> > k^T), i.e., for 
transitions related with phonon absorption, the energy- 
dependent part of r(aj) strongly decreases for increasing 
u>; for u < 0, corresponding to phonon emission, T(lu) ap- 
proaches the temperature-independent power-law depen- 
dence uj 3 . The largest r(w)'s are attained for these latter 
processes in connection with low-energy spin states. The 
contribution from the spin operators, cf. Eq.([B2]), on the 
other hand, is larger for spin states closest to the top of 
the barrier and tends to balance the changes in T(ui) and 



reduce the variations in £ 



for different states 



and for varying H z . The typical energy scale arising from 
the spin-phonon rates is found to be 10~ 5 — 10 _4 K. 



APPENDIX C: CHOICE OF BASIS 



In sections [II A and III B , we decidedly formulated the 
more general equations independent of the chosen basis 
for Tig ■ In this appendix, we consider the eigenbases of 
H.S or the (i-basis, in more detail. 

The (strong) tunneling poses a problem for the dia- 
grammatic formulation in the m-basis, but this can be 
easily solved by first diagonalizing Tis and then express- 
ing all the equations in its eigenbasis. In this d-basis, the 
kinetic equation for the diagonal and off-diagonal density 
matrix elements reads 



p{t)d,d = zZ^dd^xd^P^di^ ( C1 ) 
di 



and 



p(t)d,d> = -i{E d -E d/ )p(t) d ^ + ^2j: dd ,, dldli p(t) dudli , (C2) 



respectively. From the knowledge of (the full) p(t) we 
can again obtain, e.g., the magnetization 



M(t) =g[x B J2^2(d'\m)m{m\cC)pd,d>(t) (C3) 

d.d' m 

= gfi B ^m d ^p d ^(t) (C4) 



d.d' 



where m d d i is defined in this way as the matrix element 
of S z in the d-basis. 

When the tunneling rates dominate the kinetic equa- 
tions, the main features of the eigenstates can be under- 
stood in even simpler terms as follows. When a given 
state is off- resonant, there is essentially a one-to-one cor- 
respondence between each of the m- and d-states, i.e., 
also the d-states are localized on one or the other side of 
the barrier. Close to a resonance, two states mi and m r 
on different sides of the barrier, see Fig. ^2|, get coupled 
and form an approximate two-state system described by 



n 2 



E m , A 
A* E„ 



(C5) 



The nondiagonal elements denote the tunnel splitting as 
obtained from the diagonalization of the full spin Hamil- 
tonian; the subscripts stand for the left and right sides 
of the barrier. The eigensolutions to this are the sym- 
metric and antisymmetric combinations of the respective 
m-states 



\d s ) = a\mi) + P\m r ) 
\d a ) = (3\mi) - a\m r ) 



(C6) 
(C7) 



that extend through the barrier, see Fig. E2i The factors 
a and (3 are the normalized constants 



A 



i 

W 2 + i A i 2 



(C8) 
(C9) 



with e = \m - E r ) - y/(E t - E r ) 2 + A\A\ 2 ]. 

The biggest simplification is attained when we argue 
that, for the most values of H z , we can restrict our con- 
siderations to the diagonal elements of the density ma- 
trix. A naive justification for this concerns the stationary 
values of the density matrix elements [obtained by requir- 
ing p(t) dtd r = 0]. This leads to the immediate conclusion 
that all the off-diagonal elements between nonresonant 
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states are negligibly small. Furthermore, the nondiago- 
nal elements are also very small for any pair of resonant 
states as long as the tunnel splitting of that particular 
resonance is larger than the spin-phonon rates coupling 
these states to others, see the end of App. [b]. 

We also investigated the temporal behaviour of the off- 
diagonal elements in terms of the reduced model shown in 
Fig. ^ and the results lend support to the above conclu- 
sions. The idea of this simulation was to prepare the sys- 
tem into the state di at the initial time to, let the system 
then evolve in time according to the kinetic equation, and 
see how the off-diagonal elements p{t)d s ,d a and p{t)d a ,d 3 
behave. The reson ant pair of sta tes in the figure is similar 
to the one in Eqs. ( |C6| ) and (07) and it is coupled to two 
lower, nonresonant states di and df. The rates depicted 
in the figure are E u = E U)ii and E d w + S// irr )/2. 

The magnitudes of these rates - as compared to the tun- 
nel splitting |2A| - determine two regimes. If 2| A| 3> Ed, 
the amplitudes of the nondiagonal elements are found to 
quickly reach their maxima oc E U /|2A| and their values 
orbit around and "decay" towards the respective com- 
plex stationary values. On the other hand, according 
to the detailed-balance relation, the stationary values of 
the diagonal elements arc proportional to E u /Ed- Hence 
p{t)d,d> / p{t) d {,) oc Ed/|2A| and we can neglect the non- 
diagonal elements, if Ed -C 2|A|. In this case the ki- 
netic equation, Eqs. (CI) and (C2), becomes very simple: 
Eq. (C2) can be neglected and the rate E acquires the 
form, cf. App. |b], 



^d'd',dd ~ ±Gdd>,d>d 



nA 2 



i2pc 5 r 



AE 3 

3 f3AE _ i ■ 



(C10) 



where 



G 



dd' ,d' d / J / J rn. 4 m 3 , m 2 in, 1 
mi ,rri2 7713,7714 



(Cll) 

(d'\m4) (m.3 |g?) (d\m2) (mi \d'} 



corresponds to the vertices and contains the spin-phonon 
coupling constants, see App. [a|. In the opposite case, 
2|A| <C Ed, the nondiagonal elements do not per- 
form orbiting motion in the complex plane but increase 
motonously to roughly one half of p^ t) . In this case, the 
off-diagonal elements clearly cannot be neglected. 
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FIG. 22. Schematic of the situation discussed in the text. 
If the spin-phonon coupling denoted by the arrows and E's is 
much weaker than the tunnel coupling A between the states 
\mi) and |m r ), these states can be thought of as an effective 
two-state system. The eigenstates of the two-state system are 
depicted as dashed lines separated in energy by 2|A|. The 
subscripts of the E's correspond to those used in the text. 



For Mni2 we can attain the whole range of cases: for 
the most strongly coupled level(s) 2|A| 3> Ed, while for 
the lower levels, 2|A| <C Ed. In the former case, the di- 
agonal elements Pd(t) are sufficient in describing the sys- 
tem whereas, in the latter case, we either have to include 
the nondiagonal states or restrict our considerations to 
magnetic fields for which AE 3> Ed for all the levels, cf. 
Ref. [l(| In the text, we neglect the nondiagonal elements 
in the calculation of xi 10 ) an d m some of the analytical 
considerations but compare the two cases in section IV. 



APPENDIX D: LAPLACE TRANSFORMATION 

In this appendix, we consider the Laplace tranforma- 
tion 



f(z) = / dte-^fit) 
Jo 



(DI) 



of the kinetic equation, Eq. ([L4|). We also give another 
proof of the applicability of the Markov approximation 
in calculating the relaxation rates. 

The kinetic equation is readily transformed into 

- izp(z) - p(t = 0) = -iL p{z) + E(z)p(z) (D2) 
p(t = 0) 



P{z) 



—iz + iLo — E(z) 



(D3) 



The poles of Eq. (D3), i.e., the solutions of — izi + iL — 
E(zj) = 0, yield the exact eigenvalues to the kinetic equa- 
tion: Zi = uji + z/t;. For the slowest mode of the time 
evolution, one can consider the expansion 



E(^) «E(0)+*i 



3E(z) 



dz 



(D4) 



The prefactor of the z\ may be evaluated to be propor- 
tional to [ti ■ min{A_B, k^T, D}]^ 1 , i.e., to the maximal 
ratio between the over-barrier relaxation rate 1/ri and 
the other characteristic energy scales in the problem: 
level spacing and/or splitting AE, temperature fceT, and 
the cutoff of the phonon spectrum D (the cutoff D is in - 
troduced in order to assure that the real part of Eq. (All) 
is convergent also when we consider the nondiagonal den- 
sity matrix elements). It turns out that the actual relax- 
ation rates are several orders of magnitude smaller than 
any other energy scale and it becomes safe to approxi- 
mate 



P{z) 



pit - 0) 



= P(t = 0) 

-iz + iL - E(0) -iz-W 



(D5) 
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where S(0) has been identified as the constant £ of the 
Markov approximation above and W is defined accord- 
ingly, cf. Eq. ([15]). This approximation is valid for the 
relaxation mode and time ti 



pW(t = 0) 
— £2 — 1/ri 



(D6) 



but the Markov approximation may give erron eous re- 
sults for the faster eigenmodes for which Eq. (D5) no 
longer holds true. 



APPENDIX E: LORENTZIAN PEAK SHAPES 

The series of peaks found in the relaxation rates/times, 
cf. Fig. ^[ may be understood in terms of different relax- 
ation paths, each path with a possible tunneling channel 
giving rise to a peak - see |l8| for nice illustrations of 
the paths. In this appendix, we sketch a derivation that 
aims to show that the Lorentzian peak shapes are actu- 
ally something to be expected. 

When the spin system is somehow disturbed away 
from equilibrium and then let relax, it quickly acquires 
a metastable state, a thermal equilibrium separately on 
each side of the barrier. This initial thermalization into 
the metastable state is driven by the spin-phonon inter- 
action that can change the spin states m by ±1 or ±2. 
At a much longer time scale, the system relaxes over the 
barrier towards the real equilibrium configuration. For 
the relaxation to take place, the crucial step is the fi- 
nal transition that transfers the spin onto the other side 
of the barrier. We can distinguish two regimes in terms 
of how this critical transition takes place. In absence 
of tunneling, e.g., in off-resonance conditions, the relax- 
ation is only possible over the top of the barrier, while for 
relatively strong tunnel splitting and for resonant con- 
ditions, the dominant path is via tunneling across the 
barrier well below its top. When the tunneling is weak 
compared with the spin-phonon interaction, the tunnel- 
ing rate is the bottle neck for the relaxation to take place. 
In Mni2, this is the case for tunneling between the low- 
lying states with \m\ > 4 (for H z ps OT). However, for 
the experimentally relevant resonances, the tunneling is 
strong and takes place between the higher states. In this 
case, the spin actually oscillates back and forth through 
the barrier until it relaxes to some lower state on either 
side of the barrier. This is the case of interest here. 

In the strong-tunneling regime, the system is best de- 
scribed in terms of the (i-basis where it suffices to con- 
sider the diagonal elements of the density matrix. In 
order to get a more intuitive picture of the relaxation, 
let us consider a situation where the system has been 
prepared onto one side of the barrier and has reached 
the metastable thermal equilibrium there. This initial 
condition is convenient for two purposes: first, the relax- 
ation only proceeds into one direction and, second, the 
phonon-induced transitions on this one side of the barrier 



are accounted for by the thermal probabilities p d (tilde 
denotes the metastable state and we write just one in- 
dex for the diagonal matrix elements). Let us further 
consider relaxation via a single tunneling resonance and 
take into account the states and transition processes illus- 
trated in Fig. |2^. The system starts in the initial state di 
(localized onto the left side of the barrier, ~ rrii), is then 
activated onto resonance into either the symmetric or 
antisymmetric state, denoted by d s and d a , respectively, 
and at some point is transfered down to the final state df 
(localized onto the right side of the barrier, ~ mf). The 
subsequent intravalley relaxation is so fast that after the 
transition to df the relaxation can be considered com- 
plete. The states d s and d a extend through the barrier 
and this is the key point of the present discussion: the 
spin is transferred through the barrier in a single step, the 
rate being determined by thermal activation but also by 
the m agn etic-fie ld d ependent amplitudes a and (3, from 
Eqs. (C6) and (C7), for the two extended states to be 
on either side of the barrier. These amplitudes deter- 
mine the relative probabilities for the activation process 
to couple to the resonant states. 
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FIG. 23. Illustration of the parameters discussed in the 
text: states in the m and rf-bases, transition rates E mim / , and 



factors ±a and j3 from Eq. (C5) 



The above discussion can be formulated in the lan- 
guage of a master equation: 

p dt (t) = w Z dt ,dJdM ~ (Pd f ,d. + Z dz , ds )p ds (t) (El) 
ps \(d s \mi)\ 2 T, mitmi p d% (t) 

+ \(d s \mi)\ ,m>i 



mi .rrii Pdi \t) 
2 



PdM = ~ E 4,AP*C0 - (E d/ ,rf Q + ^ di ,da)pd a (t) (E2) 

ps \(d a \mi)\ 2 'S mil7 n i pd i (t) 
-(\(da\m r )\ 2 'E mf>mr + \ (d a \mi)\ 2 T,„ Himi )p da (t) 



\P\ 2 Y, miiTni p di (t) 



)pdM 



The approximate equalities are just a reminder that we 
have neglected, e.g., the contributions from states above 
the resonance as well as the return possibility from state 
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df. These equations can be simplified by the assumption 
^mi.mi ~ ^mf,m r which is reasonable for a pair of reso- 
nant states - as a result, the S's can be taken out of the 
brackets in the last forms of the above formulas. By fur- 
ther noting that due to normalization \a\ 2 + |/3| 2 = 1 and 
that the resulting probabilities are time independent, we 
obtain 



Pd s 
Pd a 



y 

I ^ Pdi 



,|2 ^m; ,rj . 



The ratio of the E's is just the thermal factor 
exp[— f3(Ei — Ei)], cf. detailed balance, and p~di is the 
thermal probability to be in a state with energy Ei over 
the lowest energy E- W on left hand side (for H z > 0). 
Together these yield a factor c-exp[— f3(Ei— E-io)], where 
c is a normalization constant equal to p-w which is close 
to unity for the temperatures of interest. 

In the next and final step, the relaxation rate is ob- 
tained from the knowledge of these probabilities and the 
rates to be dragged down on the right hand side of the 
barrier, i.e., 

T 1 ~ ^d f ,d e ' Pd 3 + ^d f ,d a ' Pd a 



2|a| 2 |/3| 2 £ 



rnf,m r fid s H~ \®\ t ,m r Pd a 

,„ t , mr c-e-^- E -™\ (E3) 



The exponential factor is just the effective Arrhenius fac- 
tor seen in experiments, i c£ m , mr — t^ 1 , and 



2| fl|2 



4\a\ 2 \[3 



(2|A|) 5 



(E mi -E mr )2 + ( 2 |A|)2- 



(E4) 



Here 2|A| is the tunnel splitting. It is more or less in- 
dependent of H z but £ = E mi — E mr can be tuned with 
the magnetic field. In terms of the field, the width of the 
resonant peak at its half maximum is 



SH Z = 



4IAI 



9PB\mi 



(E5) 



This sketch of a derivation introduces all the factors 
seen in experiments: the Arrhenius law with a reason- 
able prefactor t^ 1 , that depends weakly on temperature 
and the particular resonance, see the two last paragraphs 
of App. and peaks of accelerated relaxation superim- 
posed on it. The peak heights or the relaxation rates on 
resonance are found to correspond to the Boltzmann or 
Arrhenius factor with the energy corresponding to the 
effective barrier height. The peak shape is Lorentzian as 
observed in experiment with widths given by precisely 
the tunnel splittings. 



1 J.R. Friedman et al, Phys. Rev. Lett. 76, 3830 (1996). 

2 J.M. Hernandez et al, Europhys. Lett. 35, 301 (1996). 

3 L. Thomas et al, Nature (London) 383, 145 (1996). 

4 J.M. Hernandez et al, Phys. Rev. B 55, 5858 (1997). 

5 F. Luis et al, Phys. Rev. B 55, 11 448 (1997). 

6 A. Ganeschi, D. Gatteschi, and R. Sessoli, J. Am. Chem. 
Soc. 113, 5873 (1991); R. Sessoli et al, Nature (London) 
365, 141 (1993). 

7 R. Sessoli, Mol. Cryst. Liq. Cryst. 274, 145 (1995). 

8 Some authors do report possible observations of dipolar 
effects, especially at low temperature, see e.g. Refs. [l^ and 

The division between thermally-activated and tunneling- 
enhanced relaxation depends on typical experimental time 
scales. In relaxation measurements, the magnetization first 
relaxes faster than exponentially and only asymptotically 
reaches an exponential regime. For T > 3K, almost all 
the spins relax before reaching this stage. In susceptibility 
experiments, on the other hand, the relevant time scale is 
the measurement frequency and this can keep up with the 
relaxation rate up to 5-6K. 

10 J. Villain et al, Europhys. Lett. 27, 159, (1994). 

11 P. Politi et al, Phys. Rev. Lett. 75, 537 (1995). 

12 F. Hartman-Boutron, P. Politi, and J. Villain, Int. J. Mod. 
Phys. B 10, 2577 (1996). 

13 M.A. Novak and R. Sessoli, in Quantum Tunneling of Mag- 
netization, edited by L. Gunther and B. Barbara (Kluwer, 
Dordrecht, 1995). 

14 D.A. Garanin and E.M. Chudnovsky, Phys. Rev. B 56, 11 
102 (1997). 

L. Thomas, A. Caneschi, and B. Barbara, Phys. Rev. Lett. 
83, 2398 (1999). 

16 F. Luis, J. Bartolome, and J.F. Fernandez, Phys. Rev. B 
57, 505 (1998). 

17 A. Fort et al, Phys. Rev. Lett. 80, 612 (1998). 

18 M.N. Leuenberg er and D. Loss, Eu rophys. Lett. 46, 692 
(1999); preprint |co"nd-mat/9907154| . 

19 N.V. Prokof'ev and P.C.E. Stamp, Phys. Rev. Lett. 80, 
5794 (1998). 

20 T. Ohm, C. Sangregorio, and C. Paulsen, Eur. Phys. J. B 
6, 195 (1998); J. Low Temp. Phys. 113, 1141 (1998). 



D.A. Garanin, E.M . Chudnovsky, and R. Schilling, preprint 



cond-mat/9911055 

M. Al-Saqer et al, preprint |cond-mat/9909278 



' I. Tupitsyn and B. Barbara, preprint cond-mat/0002180. 
T. Pohjola and H. Schoeller, in proceedings of LT22, to 
appear in Physica B. 

' A.D. Kent it et al., preprint cond- mat/9908348 - submit- 
ted to Europhys. Lett.; Y. Zhong and M. Sarachik report 
observing side resonances also in relaxation measurements, 
private communication - to be published. 

' A.L. Barra, D. Gatteschi, and R. Sessoli, Phys. Rev. B 56, 
8192 (1997). 

Y. Zhong et al, J. Appl. Phys. 85, 5636 (1999). 

R. Sessoli et al, J. Am. Chem. Soc. 115, 1804 (1993). 



21 



This bears analogy to the related molecular magnet Fes, 
which has hard, medium, and easy axes - the interference 
effects in the tunnel splittings appear when the external 
magnetic field is oriented parallel to the hard axis. 
Similar results were recently reported in Ref. see Ref. ^f] 
for a detailed analysis of the interference effects. 
M.N. Leuenberger and D. Loss, to be published. 
In Ref. ^ the authors calculate the energy spectrum for 
Mni2 starting from the twelve atomary spins and find 
different level splittings to those predicted by the simple 
S — 10 model. The discrepancy varies between different 
pairs of states and as ratio varies between 0.8 to 5 from 
higer to lower states, respectively. Such details do not affect 
the essence of this paper and may be inserted as parameters 
into the calculation. 



W. Wernsdorfer, R. Sessoli, and D. Gatteschi, Europhys. 
Lett. 47, 254 (1999). 

D.A. Garanin, J. Phys. A 24, L61 (1991). 
H. Schoeller and G. Schon, Phys. Rev. B 50, 18436 (1994); 
J. Konig et at, Phys. Rev. B 54, 16820 (1996). 
Rigorously taken, instead of p(t') ~ p(t), we should 
use pit') = e ms(t " t ' ) p(t)e" iWs(t " t ' ) and insert this into 



Eq. (|l4[). Combining the exponential factors yields the ex- 
pressions of Refs. |4l] and |l£|. This procedure complicates 
the correspondence of the indices in the equations and in 
the diagrams, and, as this is. found to make very little dif- 
ference to the final results,tll we neglect this step here. 
We performed the calculation with both kinds of Markov 
approximations and the results did not show any significant 
differences. It is also noteworthy that when using only the 
diagonal states of the eigenbasis the two approximations 
become equal. 

We have solved the problem numerically for the full den- 
sity matrix in both bases, i.e., either first diagonalizing the 
spin part or taking the tunneling matrix elements into ac- 
count to all orders in the m-basis. The results are of course 
equivalent. 

R. Kubo, M. Toda, and N. Hashitsume, Statistical Physics 
II, Nonequilibrium Statistical Mechanics, Springer, 1985. 
E. Chudnovsky and J. Tejada, Macroscopic Quantum Tun- 
neling of the Magnetic Moment, Cambridge, 1998, ch. 6. 
K. Blum, Density Matrix Theory and Applications, 2nd edi- 
tion, Plenum Press, 1996, ch. 8. 
G. Belessa et al, Phys. Rev. Lett. 83, 416 (1999). 
The discussion on the decoherence is relevant for the low- 
lying states with small tunnel splittings. 
W. Wernsdorfer and R. Sessoli, Science 284, 133 (1999). 
D. Loss, D.P. DiVincenzo, and G. Grinstein, Phys. Rev. 
Lett. 69, 3232 (1992); J. von Delft and C.L. Henley, Phys. 
Rev. Lett. 69, 3236 (1992). 
A. Garg, Europhys. Lett. 22, 205 (1993). 
W. Wernsdorfer, contribution to the proceedings of LT22, 
Helsinki, to appear in Physica B; private communication 
with W. Wernsdorfer. 

W. Wernsdorfer et al, Phys. Rev. Lett. 82, 3903 (1999); 




22 



